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Direct numerical simulation of turbulent 
condensation in clouds 

By R. Paoli and K. Shariff f 


1. Motivation and objectives 

Clouds are responsible for precipitation, scattering and absorption of electromagnetic 
energy, and cloud particles serve as sites for chemical reactions. To properly incorporate 
these processes in cloud-resolving codes or climate models, the crucial quantity required 
is the droplet size spectrum n(r) within a computational cell, defined such that ndr 
gives the number of particles in the cell having radii between r and r + dr. The role of 
turbulence in clouds has been the object of intensive investigation in the atmospheric 
science community over the last forty years. (Consult the reviews by Pinsky & Khain 
1997, Vaillancourt & Yau 2000, and Shaw 2003). Turbulence has three effects on the 
development of the droplet size spectrum: (i) During the initial phase of condensational 
growth it causes each particle to experience a different fluctuating supersaturation as it is 
transported and this leads to a broadening of the spectrum (Cooper 1989). (ii) Turbulence 
causes particles to cluster when they become large enough that particle inertia becomes 
important; this clustering eventually influences the size spectrum during the coagulation 
phase, (iii) Finally, during the coagulation phase, turbulence increases the number of 
particle encounters. The present work aims to further elucidate the first of these effects 
with a view to developing sub-grid models for it. In the atmospheric sciences literature, 
this effect is referred to as “stochastic condensation.” 

The first step in current approaches for modeling stochastic condensation (see for 
example Khvorostyanov & Curry 1999a) is to write down the transport equation, 

df (dfuj) _d(ff) 

dt dxj dr 3 1 ' ; 

for the particle size distribution /(x,r, t). Here f(x,r,t)dxdr gives the number of par- 
ticles in the four dimensional phase space volume dxdr at (x, r). The next step follows 
the standard approach for Reynolds averaging. Decomposing variables into mean and 
fluctuating parts, f = f + f',r = r + r', and averaging the entire equation leads to the 
appearance of covariances like u'-f or r'f which are modeled using a mixing length ap- 
proach. Khvorostyanov & Curry (19996) were able to obtain a solution for / in terms of 
gamma functions which is an attractive result as many observed cloud-droplet distribu- 
tion can be fitted with gamma distributions (Shaw 2003). One of the assumptions made 
by those authors was to use a linear relation between droplet growth rate and supersat- 
uration fluctuations, f' ~ S', based on the analysis by Srivastava (1989) and heuristic 
considerations about the equivalence between Brownian motion of small particles and 
the motion of turbulent eddies. 

In this brief, we investigate the turbulent condensation of a population of droplets by 
means of direct numerical simulation. To that end, a coupled Navier-Stokes/Lagrangian 
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solver is used (see Paoli et al. 2004) where each particle is tracked and its growth by 
water vapor condensation is monitored exactly. The main goals of the study are to find 
out whether turbulence broadens the droplet size distribution, as observed in in situ 
measurements. The second issue is to understand if and for how long a correlation between 
the droplet radius and the local supersaturation exists for the purpose of modeling sub- 
grid scale microphysics in cloud- resolving codes. 

This brief is organized as follows. In Sec. 2 the governing equations are presented, 
including the droplet condensation model. The implementation of the forcing procedure 
is described in Sec. 3. The simulation results are presented in Sec. 4 together with a 
sketch of a simple stochastic model for turbulent condensation. Conclusions and the 
main outcomes of the study are given in Sec. 5. 


2. Governing equations 

The dimensionless forced compressible Navier-Stokes equations are 

(2.1) 
(2.2) 

(2.3) 

(2.4) 

Here p is the density, u = [ui,u 2 ,U 3 ] is the velocity vector, p is the pressure, E is the 
total energy, Q = [Qi, Q 2 , Q 3 ] is the heat flux vector given by Fourier’s law, C p being the 
non-dimensional specific heat at constant pressure. The quantity is the shear stress 
tensor, and Y n ,n = 1,2 are scalar fields, Y± being the temperature T and Y 2 the vapor 
mass fraction Y v . while c o n are physical source terms for Y n . In the case of water vapor 
(Y 2 =Y v ), lo 2 = u> v represents the removal of vapor due to condensation. Quantities are 
made non-dimensional using reference values: p ie f for the density, a re f for the velocity, p re f 
for the pressure, Z re f for length, T re f for temperature, p ie { for the dynamic viscosity, and 
Cp, r ef for the specific heat. The Reynolds number is defined as Re = o re f l ie t / (prei / Pret) ■ 
The Schmidt number is defined as Sc n = p ie {/ p ie tD n , where D n is the diffusivity of scalar 
n. For atmospheric conditions vapor diffusivity (given in Pruppacher & Klett 1997, p. 503 
is not too different from thermal diffusivity and so we chose Sc = 0.75 for both vapor and 
temperature. Finally, the terms pfj,pW, pqi , and pq 2 in (2.2)-(2.4) represent forcing 
of momentum, total energy, temperature and vapor, respectively and are described in 
detail in the next section. 

The use of compressible equations to simulate atmospheric clouds requires some com- 
ment since use of the incompressible or Boussinesq equations is more common. The rms u' 
of velocity fluctuations in clouds is of the order of 1 to a few meters per second (MacPher- 
son & Isaac 1977). For u' = 3 m/s and speed of sound a = 300 m/s, the turbulence Mach 
number M t = u' ja = 0.01. If we were to choose a similar M t for the simulations, the 
CFL criterion would lead to a time step that is very small compared with the flow time. 
Instead we chose M t ~ 0.05 for the simulations. This leaves the velocity field essentially 
incompressible, but allows one to run with an affordable time step. However, in weakly 
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compressible flow, temperature fluctuations are induced by pressure fluctuations and are 
proportional to Mf. Furthermore the mechanisms for temperature fluctuations in the 
atmosphere are not unrelated to flow induced pressure fluctuations. Thus we have no 
right to derive the temperature from the energy E which obeys (2.3). For this reason, 
temperature is obtained independently by solving a scalar transport equation for Y\ = T. 

2.1. Particle treatment 

Due to their small size (less than about ten microns during the condensation phase) , the 
relaxation time t p = 4 p v r 2 / 18p, of particles is short compared to the smallest turbulence 
time scale. They can then be treated as tracers which follow the gas according to 

^=u(x p ,f). (2.5) 

Gas variables at a particle location x p are estimated by linear interpolation, using the 
values at the nodes of the surrounding cell (see Boivin et al. 1998 for details). For the law 
of droplet growth by condensation we start with equation (15.74) in Seinfeld & Pandis 
(1997) and ignore the solute contribution to the equilibrium vapor pressure over the 
drop. We retain the latent heat term as well as the curvature contribution (Kelvin effect) 
and the kinetic corrections to effective vapor and temperature diffusivities. (We verified 
a posteriori that the last two contributions are negligible for drops > 1/zm). We are thus 
left with 

dr _ a(r,T)S . . 

dt r [ J 

where S = Y v — Y*(T) is the local supersaturation with respect to water and a(r,T) 
is a coefficient that depends on temperature and droplet size. Saturation conditions are 
obtained from the fit by Sonntag (1994) 

p s v = p X ® = exp( oi In T + 02 T 1 + 03 + 04 T + 05 T 2 ) (2-7) 

where the mass fraction Y V 3 * S (T) and molar fraction X£(T) at saturation are related by 
Y v s = X S J{X S V + (1 - X s u ) W air /W v ), with W air /W v = 28.85/18.01 = 1.6. 

2.2. Numerical method 

The gas transport equations (2.1)-(2.4) are discretized in physical space by means of 
the sixth order compact scheme by Lele (1992). They are advanced in time together 
with (2.5)— (2.6), using a 3rd order Runge-Kutta scheme. Periodic boundary conditions 
are used in the three directions of the computational cube which has sides Lb ox = 2 n. 
The code is parallelized using domain decomposition with MPI as the communication 
protocol. 


3. Turbulence forcing 

In order to obtain a statistically stationary velocity field, a body force f(x,f) = 

[/i(x, t), / 2 (x, t), / 3 (x, t)] is applied to the momentum equations. Physically, this rep- 
resents the effect of processes such as turbulence production by shear and the turbulent 
cascade which occur on scales larger than the computational box. To preserve univer- 
sality of the smallest scales of turbulence, only the (low) vector-valued wavenumbers k 
within the sphere k = |k| < kj are forced. The force is represented as a finite Fourier 



308 R. Paoli AND K. Shari ff 

series (see Eswaran & Pope 1988a) 

f(x,i)= J2 e ik ' x f(M) (3.1) 

k<k f 


where f(k, t) = [A (k, f), \f 2 k, f), / 3 (k, f)] and /j(k, t) = ff (k, t) +i fj (k, t) is the Fourier 
transform corresponding to k = ko[l,m,n\ (fco = 2-7T / Lb ox = 1 being the fundamental 
wavenumber), while the summation is intended over the three components, l = —N/2 + 
1, N/2; to = —N/2 + 1, N/2; n = —N/2 + 1, N/2, with the condition k = ko(l 2 + 
to 2 + n 2 ) 1 / 2 < kf. At each wavenumber, f(k, t) is obtained from the divergence-free 
projection 

f (k, t) = g(k, t) — (k ■ g(k, t)/k 2 ) k, (3.2) 

of the three-dimensional stochastic process g(k,f) = [gi (k, t) , 52 (k, £),g3(k, £)] which is 
composed of six independent Uhlenbeck-Ornstein (UO) processes corresponding to the 
real and imaginary parts of gj (k, t ) = gf(k,t) + i gj(k,t), j = 1 , . . . , 3. These processes 
are governed by the Langevin stochastic differential equation (Lemons 2002), 


dgf J { k, t) = gf J { k, t + dt) - gf J { k, t) 



dt + N/ +dt (0,l) 



(3.3) 


where 77 and af are, respectively, the autocorrelation time and standard deviation of all 
processes, while N* +dt (0, 1) is a normally distributed random number with zero mean 
and unit variance associated with the time interval [t;t + dt\. From (3.3) it can be easily 
shown that each UO process has zero ensemble- mean and is exponentially auto-correlated 
with time lag 77, i.e. 

(gf J (M))=0 (3.4) 

(sf'Msf \k,t + a)) = (3.5) 

Using the condition that fj be real (i.e. ff(— k) = ff( k), fj(— k) = —fj{ k)) and 
the requirement f v fj dV = 0 that the force not change the net linear momentum (i.e. 
ff( 0) = fj (0) = 0), one has for the forcing terms fj and W in (2.2) and (2.3) 


N/2 N/2 N/2 

fjfx, t) = Uf ( k > cos ( k • x ) - fj ( k > sin ( k • x )] 1 k < k f (3-6) 

l—l m—1 n—1 

3 r 3 

W(x,t) = ^/j( x ,£)uj( x ,£) - / y,£)wj(y,£)dU(y). (3.7) 

j= 1 Jv 3 = 1 

The quantity W represents the work done by the force and must be added in compressible 
flows to maintain consistency between the momentum and energy equations (see e.g. Kida 
& Orzag 1991). In the statistically stationary state the amount of work done by the forcing 
equals viscous dissipation to heat. The integral term in (3.7) removes thermal energy in 
the mean and is necessary to keep the mean internal energy in the computational box 
from monotonically increasing due to viscous heating. 


3.1. Scalar forcing 

Let us now consider the forcing terms q n in the scalar equations (2.4), and for the 
sake of clarity neglect the subscript n which simply identifies the scalar (temperature or 
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water vapor). The same stochastic forcing as described above is used in (2.4) to drive 
temperature and vapor fluctuations at the large scales. In actual clouds it is convection 
and mixing between air parcels, and radiation, etc. that drive the fluctuations. Following 
the usual approach, one can decompose q(x, t) as 

q(x,t) = ^ e lk ' x <?(k,t), with q(k,t) = q R (k,t) + iq 1 (k,t), (3.8) 

k<kf 

where q R (k,t) and q I (k, t) are chosen to be UO processes defined by 

dq R ’\k,t) = - qR ’ I{Kt) dt + N f t + dt (0, 1) J2v*~, (3.9) 

T i V T q 

Using the fact that q(x, t ) is real one obtains finally 

N/2 N/2 N/2 

q{x,t) = 2 ^2 ^2 [? J? (k,i)cos(k-x) - <? J (k, f) sin(k ■ x)] , k<k f . (3.10) 

l—l m= 1 n—1 

In (3.9), we chose r q = t$ assuming that the same process, namely large eddies, drive 
scalar and velocity fluctuations. The ampltitude parameter a q could be obtained by 
requiring that the standard deviation of the resulting scalar field (temperature or vapor 
concentration) be in the range of atmospheric values. Strictly speaking this can only be 
verified a posteriori, however, one can obtain a good prediction of the resulting variance 
by using the assumption of a particle system, as often done in the simulation of turbulent 
flows (see Pope 2000). In these kinds of methods, the flow is regarded as an ensemble of 
fluid particles each carrying a different value for the scalar Y which evolves according to 
a Langevin equation 

dY{t) = - Y ^ ~ ^ ® dt+ £[?"(*,*) +?*(*,*)] 

7~m , . 

k—1 

dq R ' 1 (k, t) = dt + Nt +dt (0, 1) M**- 

T q y T q 

where Nf indicates all forced modes. The first term in (3.11) is a classical exchange- with- 
the-mean mixing model, T m being the mixing time (usually taken to be the integral time 
scale r e of the flow) which is assumed to be equal to forcing time scale, r m = tj. From 
(3.11) one may then derive 

d(Y\ Nf/ 2 

-fa- = J2 *) + f )]) = 0 ( 3 - 13 ) 

k = 1 

M = - 2 °JL + 2 £ (Y[q R (k, t) + q 1 (k, t)]) . (3.14) 

at r m k=i 

Using (3.4)-(3.5) and (3.9), one finally gets to 


(3.11) 

(3.12) 


(Y) = const = (y) 0 

e _2t / re - - e “ 2 ^ 
^ \ T’e 



-2 t/r e 


(3.15) 

(3.16) 
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N Re u = ^2/3 K L e /L box Re\ k m axr] o Yl l (Tl) <jy 2 / (Y 2 ) 

64 2000 0.05 0.162 49 1.15 0.001 0.001 

128 3333 0.05 0.162 69 1.68 0.001 0.001 

Table 1. Flow-field parameters 




Figure 1 . (a) Evolution of the standard deviation of the three componets of the velocity field. 
Particles are introduced at t/r e = 10 when the velocity field is statistically stationary, (b) Energy 
spectra normalized by the Kolmogorov length r) and velocity u v . Solid line: 128 3 simulation; 
dashed line: 64 3 simulation; symbols: Re \ ~ 69 grid-turbulence experiments of Comte-Bellot & 
Corrsin (1971). 


which shows that ay achieves a stationary value ay^ which satisfies 


a 


q ~ 




(3.17) 


This relation can be used to design a q , knowing some characteristic values for ay^ in the 
atmosphere. For temperature (Y -» Y± = T) and vapor concentration (Y — > Y 2 = Y v ), 
Kulmala et al. (1997) report m 0.3 K and a-y Voo m 1.4 x 10 — 5 , which give, for the 
present simulations, a qi = 2.2 x 10 — 3 , and a q2 = 1.0 x 10 -6 . We verified a posteriori that 
these choices did indeed lead to the desired fluctuation levels (see next section). 


4. Simulation Results 

For each simulation a statistically stationary turbulent flow is first generated according 
to the method described above. The radius of the sphere of forced wavenumbers is kj = 
V8, giving a total of Nj = 92 forced modes. The standard deviation and autocorrelation 
time scale of the stochastic forcing were chosen to be, respectively, af = 1.339 x 10 -4 and 
Tf = 20.8 (the same as in Eswaran & Pope 1998b). The resulting turbulence statistics of 
the flow are provided in Table 1, where k max = nN/Lbox, K = 1/2 (< 7 ^ +cr^ 2 +cr^ 3 ) is the 
turbulent kinetic energy, L e = Ln is the longitudinal integral length scale, Re\ = u'Xjv 
is the Reynolds number based on the Taylor microscale A (see Pope 2000 for details). 
These statitsics were obtained by averaging over samples collected between t = 10 r e , 
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Figure 2. Evolution of the statistics of droplet radii r and supersaturation S. (a) Mean values 
(r) and (5). The three rising curves show (r) while the three falling curves show (5). (b) Standard 
deviations a r and as. The three rising curves show a r while the three fluctuating curves show 

as- Line types: , vapor equation is forced throughout (case 1); , vapor forcing is 

removed at particle insertion (case 2); . vapor equation is unforced throughout and the 

vapor field when particles are inserted is set to be random Gaussian (case 3). 


when the velocity fluctuations reach their statistically stationary value, up to about 
t = 50 r e ; see Fig. la. This figure shows that velocity fluctuations are isotropic; this 
was also confirmed by examining the two-point longitudinal and transverse correlation 
functions (not shown). The bulk of the simulations were performed on a 64 3 grid; the 
rest employed a 128 3 grid which gives a higher Re \ and allows validation of the method 
against the grid-turbulence experiments of Comte-Bellot & Corrsin (1971) (see Fig. lb). 

Scalar equations are also forced with the same method, with the means set to ( T ) = 
29277 and (Y v ) = 0.01415, respectively. The resulting standard deviations are ax — 
0.29 K and ay, ~ 1.4 x 10 -5 , which are in the range of atmospheric values reported in 
the literature (see Kulmala et al. 1997). 

At t/r e = 50, N p = 96 3 droplets are randomly distributed in the computational domain 
and time is then reset to zero. Their number was chosen to satisfy two major constraints: 
first, it must be high enough that accurate Lagrangian statistics can be obtained, and 
second, the particle spacing A = ( Vol /Np) 1 / 3 must be of the same order of the Kolmogrov 
scale as in clouds (see Sec. 1). This was indeed verified in the present DNS where A = 
2-7t/96 = 0.67 A or A = 1.8?y (for the 64 3 simulations). 

The reference length is taken to be Z re f = 2.5 cm which implies a dimensional box 
size of 15.7 cm (for a kinematic viscosty of air i/ re f = yu re f/yUref = 1-5 x 10 _5 m 2 /s, this 
also sets a reference velocity o re f = i^ef Re/l ie { = 1.2 m/s, and turbulent fluctuations 
°uj a ref — 6 cm/s). The droplets are initially monodisperse with radius ro = 5 yum. 

Figure 2 shows the time evolution of the statistics of droplet size and of the local 
supersaturation. Three cases were run with different types of forcing of the vapor field. 
Each case is shown using a different line type and is described in the caption. As the mean 
radius (r) increases by condensation, the mean supersaturation (5) decreases due to mass 
conservation. However, their standard deviations are quite different: o> increases, up to 
about 1 yum while as fluctuates about a certain value. Apart for a short transient in case 
3 which is due to the adjustment of water vapor to the turbulent flow, the insensitivity 
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Figure 3. Time evolution of droplet size distribution: (a) Case 1; (b) Case 3; , 

t = 0.25r e ; ■ f t = r*j ,t = 2r e ; ,t = 5r e . 


of the results to vapor forcing implies that in the present simulations supersaturation 
fluctuations are mainly controlled by temperature fluctuations. 

The increase in o> reflects broadening of the droplet size spectrum, as is often observed 
in cloud measurements (Shaw 2003). This is further shown in Fig. 3 where the distribution 
of particle radius is plotted at various instants. The mechanism of broadening is that each 
individual droplet absorbs a fraction of the available supersaturation, depending on the 
local thermodynamic conditions and water vapor concentration, and grows with its own 
rate, according to (2.6). Integrating the latter gives 

r 2 (t)=ro+ f 2aS(t)dt (4.1) 

Jo 

which shows that the square of the radius r depends on the time history of supersaturation 
S. In particular, for two droplets i and j and for a short time t the difference in the 
radius goes like rf — rj ~ t(Si — Sj). Equation (4.1) also indicates that r and S initially 
develop a correlation but subsequently turbulent fluctuations decorrelate them. This is 
clearly shown in Fig. 4 where the correlation coefficient C r s = {r 1 S') /a r as ( S' being 
evaluated at particle locations) is shown for the two cases with and without forcing 
the vapor transport equation. The figure shows that, at least for the present Re \ and 
dissipation rate, the correlation coefficient is > 0.5 up to t = ti/ 2 = 5 r e . It is interesting 
to compare ti/ 2 with the characteristic condensation time for a cloud, which coincides 
with the supersaturation absorption time (see Khvorostyanov & Curry 1999a), Tf = 
(4-7T (r) 0 D- 2 N p /Vol)~ l . In non-dimensional form, this is Tf = 167 ~ 8 r e , or Tf ~ 1.6ti/ 2 
for the current simulations. Figure 5 shows scatter plots of r and S at different times 
for case 1. As the previous figure, they show that S' /as ~ r' /ay up to t ~ 2 r e . This 
result may be used to model the droplet microphysics at the sub-grid scale level in cloud 
resolving LES, as mentioned e.g. by Paoli & Shariff (2003), and discussed next. 

4.1. Towards a stochastic model of condensation 
One goal of this work was to study the correlation between droplet radius and local 
supersaturation in a cloud and, eventually, to model it for application in cloud-resolving 
LES codes. The object of this section is to propose a modeling methodology. 

The concept of fluid particles has already been used in Sec. 3 to derive a simple re- 
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Figure 4. Temporal evolution of the radius-supersaturation correlation coefficient: solid line, 

case 1; dashed line, case 3. 



lation between a q and The same approach is used here for a system of physical 
droplets, as they have negligible inertia and exactly follow fluid particles. Thus, averag- 
ing over droplets is equivalent to volume averaging provided the number of droplets is 
“sufficiently” high. This has been verified for the present simulations. Assuming a statis- 
tically stationary temperature forcing and neglecting water vapor forcing, one can use a 
Lagrangian formulation to describe temperature and vapor evolution around a physical 
droplet, 



dt — Anap w n p r 2 r 




(4.2) 
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Of O r 

Figure 5. Evolution of the scatter plot of droplet radius and supersaturation for case 1 
(r' = r — (r ) , S’ = S — ( S }): (a) t = 0.25 r e ; (b) t = r e ; (c) t = 2r e ; (d) t = 5r e . 


dT = —— — dt + — N* +dt ( 0, 1) (4.3) 

where n p is the droplet number density and p w is the water density. Expanding the 
supersaturation 5 in (2.7) around the ambient temperature T x , one gets after some 
algebra 


Y S (T) ~ r^Too) + 


(T-r 00 ) = y,(r 00 )[n- V00 (r-r 00 )] (4.4) 


with (pm = ai T^ 1 — ai T^ 2 + a\ + 2a$ T x . Substituting the latter into Eqs. (4.2) and 
(4.3) and rewriting for completeness the equation for the growth of droplet radius gives 

S — (S) I 

dS = dt + <^oo n (Too) ^24^ dt — 4w ap w n p r 2 r (4.5) 

T e V 

, otS , , , 

dr = — dt (4.6) 

r 

Taking the average of these equations over a population gives 

= -4tt ap w n p (S r) = -4nap w n p (( 5 ) (r) + C Sr o r o s ) (4.7) 

dt dt dt ' ' v 1 


The problem then is to close, using information from the DNS, the covariances that 
appear when taking higher moments of r and 5 (such as, for example, the correlation 
between r' and S'). This represents the object of our current research. 


5. Conclusions 

In this brief, a direct numuerical simulation of the turbulent condensation of a popu- 
lation of droplets was perfomed using a coupled Navier-Stokes/Lagrangian code. Forcing 
of the small wavenumbers was used to sustain velocity, vapor, and temperature fluctua- 
tions. The resulting supersaturation fluctuations were responsible for the broadening of 
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the droplet size distribution in agreement with in situ measurements. Finally, a sketch 
of a possible approach for modeling the correlation between droplet size and supersatu- 
ration for use in cloud resolving LES was presented. 

Computational resources were provided by the NAS Supercomputing Division at NASA 
Ames Research Center which is gratefully aknowledged. The first author wishes to thank 
the Center for Turbulence Research for its hospitality during his post-doctoral stay at 
NASA Ames and Stanford University. 
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